Nucleation of vortex arrays in rotating anisotropic Bose-Einstein condensates 
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The nucleation of vortices and the resulting structures of vortex arrays in dilute, trapped, zero- 
temperature Bose-Einstein condensates are investigated numerically. Vortices are generated by 
rotating a three-dimensional, anisotropic harmonic atom trap. The condensate ground state is 
obtained by propagating the Gross-Pitaevskii equation in imaginary time. Vortices first appear at 
a rotation frequency significantly larger than the critical frequency for vortex stabilization. This 
is consistent with a critical velocity mechanism for vortex nucleation. At higher frequencies, the 
structures of the vortex arrays are strongly influenced by trap geometry. 
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Since the experimental achievement of Bose-Einstein 
condensation (BEC) in confined alkali gases jO-pl, the 
possibility of generating vortices in confined weakly- 
interacting dilute Bose gases has been intensively stud- 
ied |5H12|. While theoretical investigations of stability 
have generally been restricted to the case of a single vor- 
tex [0-071, the proposed experimental techniques may 
induce several vortices simultaneously Jig ]. Under ap- 
propriate stabilizing conditions, such as a continuously 
applied torque, these vortices would form an array akin 
to those obtained in rotating superfluid helium J13 Ell . 

A standard approach used to 'spin-up' superfluid he- 
lium is to rotate the container at an angular frequency 
Q. Aside from significant hysteresis effects [ p2|]23| , vor- 
tices tend to first appear at a frequency Q v , whose value 
is comparable to the critical frequency fi c at which the 
presence of vortices lowers the free energy of the inter- 
acting system ]24| . Energy minimization arguments have 
also yielded vortex arrays that are very similar to those 
observed experimentally (2^j2^]. Despite these successes, 
the mechanisms for vortex nucleation by rotation remain 
poorly understood; important factors are thought to in- 
clude presence of a normal fluid, impurities, and surface 

roughness H,|l|6). 

It has been suggested that vortices may be similarly 
generated in the dilute Bose gases by rotating the trap 
about its center p4 18|. Evidently, a harmonic potential 



can transfer angular momentum to the gas only if it is 
anisotropic in the plane of rotation. While vortices in 
such a system at zero temperature have been shown to 
become energetically stable for Q > Q c |T^JlJ,|lj, the 
particle flow could remain irrotational at these angular 
frequencies since there exists an energy barrier to vortex 
formation |14| . Suppression of this barrier could be in- 
duced by application of a perturbing potential near the 
edge of the confined gas, as has been simulated in the 
low-density limit Eq|. One of the primary motivations 
for the present work, however, is to determine if there 
exists any intrinsic mechanism for vortex nucleation in 
a dilute quantum fluid that is free of impurities, surface 
effects, and thermal atoms. We find that vortices can 



indeed be generated by rotating Bose condensates con- 
fined in an anisotropic harmonic trap. The value of f2„ at 
which vortices are spontaneously nucleated is somewhat 
larger than fl c . For Q > Q v multiple vortices appear 
simultaneously, in patterns that depend upon the geom- 
etry of the trap. 

The dynamics of a dilute Bose condensate at zero 
temperature are governed by the time-dependent Gross- 
Pitaevskii (GP) equation p7fl . Previous simulations 
of the GP equation have demonstrated that vortex- 
antivortex pairs or vortex half-rings can be generated 
by superflow around a stationary obstacle ^HgJ^ or 
through a small aperture p{|. In ||, the vortex pairs 
form when the magnitude of the superfluid velocity ex- 
ceeds a critical value which is proportional to the lo- 
cal sound velocity; recent experimental results support 
this conclusion |30| . To our knowledge, no numerical in- 
vestigation of vortex nucleation in three-dimensional in- 
homogeneous rotating superfluids has hitherto been at- 
tempted. 

The numerical calculations presented here model the 
experimental apparatus of Kozuma et al. |3l| , where 23 Na 
atoms are confined in a completely anisotropic three- 
dimensional harmonic oscillator potential. In the pres- 
ence of a constant external torque, the condensate obeys 
the time-dependent GP equation in the rotating reference 
frame: 

id T ^(r, t) = [T + Krap + Vh - nLJ V(r, r), (1) 



where the kinetic energy is T = — |V 2 , the trap po- 
tential is Vtrap = \ (x 2 + a2 y 2 + fl 2 ? 2 ), an d the Hartree 
term is Vh = 47r?7|'0| 2 . The angular momentum operator 
L z = i {yd x — xd y ) rotates the system about the z-axis 
at the trap center at a constant angular frequency fi. 
The trapping frequencies are (ui x , U) y , U) z ) = u> x (l,a,/3) 
with lj x = 2tt x 26.87 rad/s, a = \/2, and (3 = I/a/2- 
Normalizing the condensate J dr\ip(r, r)| 2 = 1 yields the 
scaling parameter n = Noa/d Xl where a = 2.75 nm is 
the s-wave scattering length for Na and N is the num- 
ber of condensate atoms. Unless explicitly written, en- 



ergy, length, and time are given throughout in scaled 
harmonic oscillator units Tiu> Xl d x = \jT%jMui x « 4.0 /im, 
and T = clTT 1 s=s 6 ms, respectively. 

The stationary ground-state solution of the GP equa- 
tion, defined as that which minimizes the value of the 
chemical potential, is found by norm-preserving imagi- 
nary time propagation (the method of steepest descents) 
using an adaptive stepsize Runge-Kutta integrator. The 
complex condensate wavefunction is expressed within 
a discrete-variable representation (DVR) |3^] based on 
Gauss-Hermitc quadrature, and is assumed to be even 
under inversion of z. The numerical techniques are de- 
scribed in greater detail elsewhere 0,p2f. The initial 
state (at zero imaginary time f = ir = 0) is taken to be 
the vortex-free Thomas-Fermi (TF) wavefunction t/jtf = 
VvMtf — Vtrapl/47r?7, which is the time-independent so- 
lution of Eq. (Q), neglecting T and L z , with chemical 
potential (j,tf = ^(15a/?7y) 2 / 5 . The GP equation for 
a given value of fi and Ao is propagated in imaginary 
time until the fluctuations in both the chemical potential 
and the norm become smaller than 10~ n . It should be 
emphasized that the equilibrium configuration is found 
not to depend on the choice of purely real initial state. 
Since the final state is unconstrained except for z-parity, 
the lowest-lying eigenfunction of the GP equation corre- 
sponds to a local minimum of the free energy functional. 

In Fig. H are depicted the condensate density, which is 
stationary in the rotating frame, as well as the conden- 
sate phase and the velocity field in the laboratory and 
rotating frames, for fi = 0.45^ and Ao — 10 6 . The den- 
sity profile at this angular frequency contains no vortices 
but is slightly extended from that of a non-rotating con- 
densate due to the centrifugal forces. The velocity field 
in the laboratory frame is given by v^ = Vy in units of 
w x d x , where tp is the condensate phase. In the rotating 
frame, v r s = v s — Clz X r. There are no closed veloc- 
ity streamlines found in Fig. 0(a). Such an irrotational 
flow V x v s = is characteristic of a superfluid, distinct 
from the related properties of vortex quantization and 
stability. The only solution of the GP equation satis- 
fying irrotational flow in a cylindrically-symmetric trap 
is v s = 0: rotating the trap is equivalent to doing noth- 
ing. The irrotational velocity field for an anisotropic trap 
is nontrivial, however. Since the density profile is inde- 
pendent of orientation, mass flow must accompany the 
rotation even though the superfluid prefers to remain at 
rest H). 

The condensate is found to remain vortex-free for an- 
gular velocities significantly larger than the expected 
critical frequency for the stability of a single vortex 



>(i) 



flc 1 ' 1l3|,p4 17|. In order to determine if irrotational con- 
figurations correspond to the global free energy minima of 
the system, vortex states are investigated by artificially 
imposing total circulation nn on the condensate wave- 
function. By winding the phase at f = by 2nn about 
the trap center, imaginary-time propagation of the GP 
equation yields the minimum energy configuration with n 



vortices if such a solution is stationary or metastable JL7) . 
The results for Nq = 10 6 and fi = 0A5u> x are summa- 
rized in Table H. At this angular frequency, states with 
n = 1,2,3 are all energetically favored over the vortex- 
free solution. The vortices in these cases are predomi- 
nantly oriented along the (z) axis of rotation, and are 
located symmetrically about the origin on the (loose) x- 
axis. The frequency chosen is too low to support the four 
vortex case fi < Q c , but is larger than the frequency 
(which may correspond to metastability) at which the 
chemical potentials for n — 3 and n = 4 cross. 




FIG. 1. The condensate density integrated down the axis 
of rotation (brightness) and the phase ip through the z — 
plane (colors in the sequence red-green-blue-red for tp — 
through 2tt) are shown for No — 10 6 and an angular frequency 
Q, = 0A5u>x applied counterclockwise. The irrotational ve- 
locity field in the laboratory (a) and rotating (b) frames of 
reference are shown as arrows, whose lengths are proportional 
to the magnitude of the velocity |v s | at the quiver tail. 

As shown in Fig. @, vortices with the same circu- 
lation k (as opposed to vortex-antivortex pairs) begin 
to penetrate the cloud above a critical angular veloc- 
ity for vortex nucleation fi„. The value of fi„ is found 
to not depend strongly on trap geometry and to de- 
crease very slowly with Ao; for Ao = 10 9 with q = 
{5,6,7}, we obtain fi„ = {0.65, 0.50, 0.36}w x ± 0.01w K , 
respectively. In contrast, the critical frequency for the 
stabilization of a single vortex in an anisotropic trap 
is approximately given by the TF expression fie s=s 
(5a/2R 2 ) ]n(R/£)uj x , where R = VVJF and £ = y/a/R 
are the dimensionless condensate radius along x and 
the healing length, respectively |IJ],[l7|]. For the pa- 
rameters considered here, the values are predicted to 
be fi c = {0.61, 0.33, 0.16}^ and are found numeri- 
cally to be n ( c 1} = {0.54, 0.29, 0.14}u> x . The number of 
vortices n v present just above £l v is found to increase 
with A ; n v = 4 and 8 for A = 10 6 and 10 7 , respec- 
tively. The value of Q u may be interpreted as the critical 
frequency fie for the stabilization of n 7J vortices. If 



= n for all Aq, then fij, ~ fi 



(n) 
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That fi^ 



decreases more slowly with Ao implies that n„ must in- 
crease with Ao. The small difference between fi c ' and 
fiy for Ao = 10 5 reflects the instability of vortex arrays 
in the low-density limit. As Aq decreases, the spacing be- 



tween successive Qc diminishes, and vanishes for Nq = 
in cylindrically-symmetric traps; for very large N , the 
spacing approaches a constant as the vortex- vortex inter- 
actions become negligible. 




FIG. 2. The condensate densities integrated down the axis 
of rotation (perpendicular to the page) are shown in the ro- 
tating frame as a function of angular velocity Q for iVo = 10 6 . 
From the top left to bottom right in raster order, are shown 
Q — 0.50^, 0.550^, 0.65^, 0.7lu x , 0.75lj x , and 0.8^. The 
density is proportional to the brightness, and the dark spots 
correspond to unit vortices. Each box is 18d^ ~ 73/im wide. 

The numerical results for f2„ suggest that the criteria 
for vortex stabilization and nucleation are different. Su- 
perflow through microapertures, or the motion of an ob- 
ject or ion through a superfluid, can give rise to vortex 
half-ring p9|]33| ] or vortex-pair |5|Jq,p|, p8|J3C| ] production 
through the accumulation of phase-slip. One might ex- 
pect similar excitations in a rotating condensate p3| , |34fl : 
vortex half-rings would be nucleated at the condensate 
surface when the local tangential velocity exceeds a criti- 
cal value. Indeed, the distinction between a half-ring and 
vortex becomes blurred in a trapped gas with curved sur- 
faces, as discussed further below. 

A crude estimate of VL U may be obtained by invok- 
ing the Landau criterion for the critical velocity v cr = 
min(u>q/q), where u q is the frequency of the mode at 
wavevector q. Such a minimum corresponds to values 
of q c at which the hydrodynamic description of the col- 



lective excitations begins to fail |13]. For a spherical 
trapped Bose gas, the crossover to a single-particle be- 
havior occurs in a boundary-layer region at the cloud 
surface whose thickness is several 5 — (2R)^ 1 ' s d x |p5| , p6| . 
Minimizing ujq/q using the dispersion relation for the pla- 
nar surface modes |37| of such a system ui 2 rs ui 2 [qR — 



dlq A (ln(gtf) - 0.15)], one obtains q c = (.R/0.3) 1 / 3 ^ 



8 1 and VL V = v cr /R 
critical frequency Q, v ~ 



R- 2 ^. Since R 



N, 
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decreases far more slowly 



than does the TF estimate for fL 



N, 
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The number- 



dependence of f2„ is in reasonable agreement with the 
numerical data. Real-time simulations further confirm 
that high-frequency oscillations of the condensate are re- 
quired for vortex production at the same £l v found using 
the imaginary-time approach. 

The above analysis does not clearly identify the in- 
stability of the surface modes with the penetration of 
vortices into the cloud, however. Further insight may be 
gained by considering the free energy F of a single vortex 
in a cylindrical trap, relative to that of the vortex-free 
state, as a function of the vortex displacement p from 
the trap center |L4J. In the TF limit, F vanishes for 
p 2 = R 2 and p 2 = R 2 — (5/217) ln(i?/£), corresponding to 
the right and left roots of the free energy barrier to vor- 
tex generation, respectively. As f2 increases, the energy 
barrier at the surface narrows but remains finite. Yet, as 
discussed above, the hydrodynamic excitations begin to 
break down at a radius /i« R — S. Vortices will there- 
fore spontaneously penetrate the cloud when the angular 
frequency exceeds 



5^2 /JT 

4R 2 /3 \£ ] 



(2) 



since the barrier effectively disappears when p < p. Thus, 
the frequencies of nucleation and penetration have the 
same number-dependence and are defined by a single crit- 
ical wavelength. Once the condensate contains vortices 
at a given Q > f2„, the functional F will again include a 
barrier to vortex penetration from the surface, reflecting 
the hydrodynamic stability of the vortex state. One may 
thus envisage a succession of multiple- vortex nucleation 
events at well-defined angular frequencies. 

The stationary configurations of vortex arrays are 
shown as a function of applied rotation in Fig. @. The 
condensate density is shown integrated down the axis of 
rotation z, in order to mimic an in situ image of the cloud. 
While the vortices near the origin appear to have virtu- 
ally isotropic cores, those in the vicinity of the surface 
are generally wider and are noticeably distorted. The 
anisotropy is due in part to the divergence of the coher- 
ence length as the density decreases, but is mostly the 
result of vortex curvature. Off-center vortices are not 
fully aligned with the axis of rotation z, since they ter- 
minate at normals to the ellipsoidal condensate surface. 
Far from the origin, the vortex structure approaches that 
of a half-ring pinned to the condensate surface. 



The symmetries of the confining potential impose con- 
straints on the vortex arrays that may be produced by 
rotating anisotropic traps. Stationary configurations are 
found to always have the inversion symmetry (x, y, z) — > 
(—x, —y, z). As shown in Fig. 0, the number of vortices is 
at least four and is even for each array; real-time simula- 
tions demonstrate that vortices with the same circulation 
are nucleated in pairs at inversion-related points on the 
surface. No vortex is found at the origin, since the odd 
number of remaining vortices cannot be distributed sym- 
metrically. At low angular velocities, therefore, the array 
tends to approximate a regular tetragonal lattice. As the 
total number of vortices increases with f2, however, a 
different pattern begins to emerge. While a triangular 
array is inconsistent with the twofold trap symmetries, 
it is more efficient for close packing; this geometry is 
favored for vortices near the rotation axis of rapidly ro- 
tating vessels of superfluid helium p9|-pl[|. If vortices 
in trapped condensates could be made sufficiently nu- 
merous, they would likely form a near-regular triangular 
array but with the central vortex absent. 

In summary, the critical frequencies for the zero- 
temperature nucleation of vortices ft v in rotating 
anisotropic traps are obtained numerically, and are found 
to be larger than the vortex stability frequencies Vt c . The 
number-dependence of Q, v is consistent with a critical- 
velocity mechanism for vortex production. The struc- 
tures of vortex arrays are strongly affected by trap ge- 
ometry, but approach triangular at large densities. 
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19.517 
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TABLE I. The chemical potential (x, ground state energy 
E and average angular momentum per particle (L z ) are given 
as a function of the number of vortices n for No — 10 6 and 
applied angular frequency Q, = 0A5uJx- The n = 4 solution 
may be metastable since /14 < ^3 while S4 > E3. 



